Energy Modelling - Main Program in Jupyter Notebook#
Project: MIRACA#
License: MIT License#
Copyright (c) 2024 miracaEU Contributors#
See LICENSE file for details.#
Import Modules#
import os
# If needed DLL Directory uncomment below line
#os.add_dll_directory(os.path.join(os.environ['CONDA_PREFIX'], 'Library', 'bin'))
import shapely
import openpyxl
import xlrd
import pandapower as pp
import pandas as pd
import sys
sys.path.append('./scripts')
from __prompt_for_countries import prompt_for_valid_countries, get_country_code, transform_iso_code
from _0_filter_network import filter_network_by_countries, connect_DC_elements
from _1_load_processing import process_and_save_profiles
from _2_gis_processing import process_regions
from pandapower.plotting import pf_res_plotly, simple_plotly
from _3_electricity_demand import distribute_load_to_nuts
from _3_2_get_fuel_price import get_fuel_price, get_co2_prices
from _3_3_marginal_costs import load_costs, process_powerplants_data
from shapely.geometry import Point
from _4_mapping_gen_and_load import create_gen_or_load
from _5_check_connections import check_connections, verify_slack_connection
from _7_time_series_and_pf_calculations import (
time_series_pf_results,
adjust_hvdc_link_capacity,
time_series_pf_results_separate_exports,
)
from _8_1_plot_NUTS_regions_and_connections import (
plot_regions,
plot_bus_nuts_connections,
plot_bus_gen_connections,
)
from _8_2_plot_generation_and_consumption import plot_power_flow_with_lines
Available EU countries: ['Albania', 'Austria', 'Belgium', 'Bosnia and Herzegovina', 'Bulgaria', 'Croatia', 'Czech Republic', 'Denmark', 'Estonia', 'Finland', 'France', 'Germany', 'Greece', 'Hungary', 'Ireland', 'Italy', 'Latvia', 'Lithuania', 'Luxembourg', 'Montenegro', 'Netherlands', 'North Macedonia', 'Norway', 'Poland', 'Portugal', 'Romania', 'Serbia', 'Slovakia', 'Slovenia', 'Spain', 'Sweden', 'Switzerland', 'Ukraine', 'United Kingdom', 'Kosovo', 'Moldova']
Prompt for countries#
# Main script logic
if __name__ == "__main__":
countries_to_filter_full, countries_to_filter = prompt_for_valid_countries()
countries_to_filter_iso = get_country_code(countries_to_filter, iso=True)
countries_to_filter = [transform_iso_code(code) for code in countries_to_filter_iso]
print()
print(" Countries to filter:", countries_to_filter)
print()
print(" Countries to filter (ISO codes):", countries_to_filter_iso)
start_time = "2015-01-01 00:00:00"
end_time = "2016-01-01 01:00:00"
Countries to filter: ['SI']
Countries to filter (ISO codes): ['SI']
Import Pandapower Network Structure#
net = pp.create_empty_network()
# Define the folder path
folder = "Data\\Net_structure_data"
# Import each DataFrame from CSV files using the correct folder path
net.bus = pd.read_csv(os.path.join(folder, "bus_data.csv"))
net.bus_geodata = pd.read_csv(os.path.join(folder, "bus_geodata.csv"))
net.trafo = pd.read_csv(os.path.join(folder, "trafo_data.csv"))
net.line = pd.read_csv(os.path.join(folder, "line_data.csv"))
net.line_geodata = pd.read_csv(os.path.join(folder, "line_geodata.csv"))
print()
print("Import completed!")
print()
net.line.drop(columns=["std_type"], inplace=True)
net.trafo.drop(columns=["std_type"], inplace=True)
net.line.drop(columns=["type"], inplace=True)
# Filter only specific country/countries
from _0_filter_network import filter_network_by_countries, connect_DC_elements
# Filter the network for the specified countries
net = filter_network_by_countries(net, countries_to_filter_iso)
print()
print("Filtered network successfully!")
print()
Import completed!
Filtered network successfully!
1 Importing processing data (load profile)#
# Load the dataset
data = pd.read_csv(os.path.join(os.getcwd(), 'Data\\time_series_60min_singleindex.csv'), index_col=0, parse_dates=[0], date_format="%Y-%m-%dT%H:%M:%SZ")
# Get the first and last utc_timestamp from the index
first_timestamp = data.index[0] # First value in the index
last_timestamp = data.index[-1] # Last value in the index
print()
print("First UTC Timestamp:", first_timestamp)
print("Last UTC Timestamp:", last_timestamp)
print()
# Define default start and end times
default_start_time = "2015-01-01 00:00:00"
default_end_time = "2015-02-01 01:00:00"
# Prompt the user for start and end times, using defaults if input is empty
start_time = input(f"Enter the start time (default: '{default_start_time}'): ") or default_start_time
end_time = input(f"Enter the end time (default: '{default_end_time}'): ") or default_end_time
print()
# Print the selected times
print("Start time:", start_time)
print("End time:", end_time)
print()
load_profiles = process_and_save_profiles(data, countries_to_filter_iso, start_time=start_time, end_time=end_time)
First UTC Timestamp: 2014-12-31 23:00:00
Last UTC Timestamp: 2020-09-30 23:00:00
Start time: 2015-01-01 00:00:00
End time: 2015-02-01 01:00:00
2 Importing NUTS3 Data - Proccesing Regions#
# Call the function for 2019 data
processed_regions = process_regions(
nuts_path=os.path.join(os.getcwd(), 'Data\\NUTS_RG_01M_2024_4326.shp'),
uk_path=os.path.join(os.getcwd(), 'Data\\UK_NUTS_2021.gpkg'),
gdp_data_path=os.path.join(os.getcwd(),'Data\\estat_nama_10r_3gdp$defaultview_filtered_en.csv'),
pop_data_path=os.path.join(os.getcwd(),'Data\\estat_nama_10r_3popgdp$defaultview_filtered_en.csv'),
country_list=countries_to_filter,
period=2019,
nuts_level="NUTS3"
)
3 Households and Industrial Demand#
# Call the function
distributed_load = distribute_load_to_nuts(
load_df=load_profiles,
regions=processed_regions,
distribution_key={"gdp": 0.6, "pop": 0.4}
)
print(distributed_load[0:2])
id LEVL_CODE country name NUTS_NAME MOUNT_TYPE URBN_TYPE \
1318 SI031 3.0 SI Pomurska Pomurska NaN 3.0
1319 SI032 3.0 SI Podravska Podravska NaN 3.0
COAST_TYPE geometry level1 \
1318 NaN POLYGON ((16.24176 46.87629, 16.26450 46.87200... SI0
1319 NaN POLYGON ((15.78642 46.70747, 15.79946 46.69724... SI0
level2 level3 gdp pop load_factor load
1318 SI03 SI031 15700.0 114.38 0.061250 0.0
1319 SI03 SI032 18700.0 324.83 0.109065 0.0
3.1 Retrieve and build monthly prices#
# Retrieves monthly fuel prices
fuel_price = get_fuel_price("Data\\fuel_prices.xlsx")
3.2 Adding Costs#
# Costs from: https://github.com/PyPSA/technology-data/blob/master/outputs/costs_2020.csv
costs = pd.read_csv(os.path.join(os.getcwd(), 'Data\\costs_2020.csv'))
path = os.path.join(os.getcwd(), 'Data\\costs_2020.csv')
costs = load_costs(path, 1)
3.3 Dynamic prices (fuel, co2)#
start_interval = pd.to_datetime(start_time)
end_interval = pd.to_datetime(end_time)
# Filter fuel prices within start and end time
filtered_prices = fuel_price.loc[(fuel_price.index >= start_interval) & (fuel_price.index <= end_interval)]
# Check if we have more than two rows
if len(filtered_prices) > 2:
filtered_prices = filtered_prices.iloc[:-1] # Exclude the last row
monthly_fuel_price = filtered_prices.mean().to_frame().T # Converts mean result into a row format
else:
monthly_fuel_price = filtered_prices.iloc[[0]] # If two or fewer rows, take the first row
path = os.path.join(os.getcwd(), 'Data\\powerplants.csv')
plants = process_powerplants_data(path, costs, monthly_fuel_price)
data_folder = os.path.join('Data\\emission-spot-primary-market-auction-report-2012-2024-data')
co2_prices = get_co2_prices(data_folder, start_interval, end_interval)
co2_prices.index = pd.to_datetime(co2_prices.index)
# Calculate the monthly average price
monthly_co2_prices = co2_prices.resample("ME").mean()
plants["adjusted_marginal_price"] = plants["CO2 intensity"] * monthly_co2_prices.iloc[0] + plants["marginal_cost"]
2015
4 Mapping Locations of Supply and Demand#
plants["country_code"] = plants["country"].apply(get_country_code)
plants = plants[plants['country_code'].isin(countries_to_filter)]
# Create a 'geometry' column from lat and lon data
plants['geometry'] = plants.apply(lambda row: Point(row['lon'], row['lat']), axis=1)
# Add a new column 'zone' by applying the get_country_code() function to the 'country' column
plants["zone"] = plants["country"].apply(lambda x: get_country_code(x, True))
print()
user_input_gen = input("Do you want to distribute plants to their closest bus (without voltage optimization)? (yes/no): ").strip().lower()
print()
distribute_gen = True if user_input_gen in ["yes", "true"] else False
if distribute_gen:
print("Plants will be distributed to the closest bus without voltage optimization.")
else:
print("Using voltage optimization and real substations for generator mapping.")
print()
user_input_load = input("Do you want to distribute loads evenly within NUTS regions (without voltage optimization)? (yes/no): ").strip().lower()
print()
distribute_load = True if user_input_load in ["yes", "true"] else False
if distribute_load:
print("Consumers will be distributed between the NUTS regions without voltage optimization.")
else:
print("Using voltage optimization and closest real substations for load mapping.")
Plants will be distributed to the closest bus without voltage optimization.
Consumers will be distributed between the NUTS regions without voltage optimization.
4.1 Mapping powerplants#
# Call function to create gen
create_gen_or_load(gdf=plants, Type="gen", net=net, distribute_gen=distribute_gen)
power_sum = net.gen["p_mw"].sum()
load_bus_len = len(net.bus)-len(net.gen)
power_per_load = power_sum/load_bus_len
print()
print("Total power of generators: ", power_sum, " MW", "(",power_sum/1000, " GW)")
print()
Total power of generators: 3844.7673661369117 MW ( 3.844767366136912 GW)
4.2 Mapping Loads Based on NUTS 3 Regions - p.u.#
loads_NUTS3 = distributed_load
loads_NUTS3["zone"] = loads_NUTS3["country"].apply(lambda x: get_country_code(x, True))
# Call function to create loads based on NUTS-3 regions
create_gen_or_load(gdf=loads_NUTS3, Type="load", net=net, distribute_load=distribute_load)
5 Assign Slack Bus and Check Connections#
# Dynamically assign the slack bus to the generator with the highest max_p_mw
biggest_unit_idx = net.gen["p_mw"].idxmax() # Get the index of the generator with the highest max_p_mw
net.gen["slack"] = False # Set all generators' slack status to False
net.gen.loc[biggest_unit_idx, "slack"] = True # Assign slack bus to the biggest unit
5.1 Save Net Elements and Check Connections - Open Buses#
# Save net elements:
net_bus = net.bus
net_load = net.load
net_gen = net.gen
net_line = net.line
net_trafo = net.trafo
result = check_connections(net_bus, net_load, net_gen, net_line, net_trafo)
# Find the index of the current slack generator
slack_gen_index = net.gen.loc[net.gen["slack"] == True].index[0]
verify_slack_connection(result, slack_gen_index, net)
6 DC OPF Settings#
net.line["max_loading_percent"] = 100
net.trafo["max_loading_percent"] = 80
for i, cost in enumerate(net.gen["cost_per_mw"]):
pp.create_poly_cost(
net,
element=i, # Generator index
et="gen", # Element type (generator)
cp1_eur_per_mw=cost # Marginal cost from `net.gen.cost_per_mw`
)
# Connect DC equivalent elements
connect_DC_elements(net)
6.1 Check Connections After DC Equivalents#
# Step 1: Identify valid and invalid lines
invalid_lines = net.line[
(~net.line["from_bus"].isin(net.bus.index)) |
(~net.line["to_bus"].isin(net.bus.index))
]
if not invalid_lines.empty:
print(f"Warning: {len(invalid_lines)} invalid lines found. Proceeding with filtering.")
# Step 2: Keep only valid lines
valid_lines = net.line[
(net.line["from_bus"].isin(net.bus.index)) &
(net.line["to_bus"].isin(net.bus.index))
]
net.line = valid_lines.reset_index(drop=True)
# Step 3: Update line_geodata to match net.line
net.line_geodata = net.line_geodata.loc[valid_lines.index].reset_index(drop=True)
print("Dropped invalid lines and updated line_geodata!")
# Find missing buses referenced in from_bus and to_bus
missing_from_bus = set(net.line["from_bus"]) - set(net.bus.index)
missing_to_bus = set(net.line["to_bus"]) - set(net.bus.index)
print("Missing from_bus indices:", missing_from_bus)
print("Missing to_bus indices:", missing_to_bus)
else:
print("No missing data found. Skipping filtering process.")
No missing data found. Skipping filtering process.
7 Time Series and Repeated Calculations#
load_profiles = load_profiles.dropna(subset=[load_profiles.columns[0]])
time_series_short = load_profiles[0:10]
def filter_load_profiles(regions, time_series):
filtered_columns = set() # Use a set to avoid duplicates
for _, region in regions.iterrows():
country_code = region["country"]
expected_prefix = "GB_UKM_load_actual" if country_code == "UK" else f"{country_code}_load_actual"
matching_columns = [col for col in time_series.columns if col.startswith(expected_prefix)]
filtered_columns.update(matching_columns) # Add unique values to the set
filtered_columns = list(filtered_columns) # Convert set back to list
filtered_time_series = time_series[filtered_columns]
return filtered_time_series
load_profiles_ = filter_load_profiles(regions=processed_regions, time_series=load_profiles)
load_profiles_short = filter_load_profiles(regions=processed_regions, time_series=time_series_short)
def is_large_system(net, threshold=1000):
if len(net.bus) >= threshold:
return load_profiles[:1] # Return shortened time series for large systems
return time_series_short # Keep the original time series for smaller systems
time_series_short = is_large_system(net)
# If you want to include the whole interval uncomment below code:
# time_series_short = load_profiles
7.1 Joined Results - Total Load/Generation and Maximum Line Loading#
# Replace NaN values in p_mw with 0
net.gen["p_mw"] = net.gen["p_mw"].fillna(0)
net.gen["min_p_mw"] = net.gen["min_p_mw"].fillna(0)
net.gen["max_p_mw"] = net.gen["max_p_mw"].fillna(0)
adjust_hvdc_link_capacity(net)
def handle_pf_results(
mode, net, time_series_short, processed_regions, pp, variable=False
):
"""Attempts to run power flow (PF) analysis for the given mode. Returns True if successful, False otherwise."""
try:
pf_results = time_series_pf_results(
mode, time_series_short, net, processed_regions, pp, True, variable
)
# Check if all values in the results are NaN (indicating a failure)
if pf_results.isna().all().all():
print(f"All values in pf_results are NaN, skipping plot for {mode}.")
return False
# Check if any line loading exceeds 100%, indicating overloading
elif (pf_results["max_line_loading_%"] > 100).any():
print(f"Line overloading detected for {mode}.")
return False
else:
print(pf_results)
pf_res_plotly(net) # Generate visualization
print(f"{mode} converged successfully.")
return True
except Exception as e:
# Catch any unexpected errors and log them
print(f"Error in {mode}: {e}")
return False
7.1.2 Normal AC powerflow calculation#
# **Step 1: Perform AC power flow analysis**
handle_pf_results("ac_pf", net, time_series_short, processed_regions, pp)
converged = False
Updating loads and generators (solar, wind) profiles - this might take a while depends on countries input size
time_step total_load_ref_MW total_load_res_MW \
0 2015-01-01 00:00:00 1045.47 1045.47
1 2015-01-01 01:00:00 1045.47 1045.47
2 2015-01-01 02:00:00 1004.79 1004.79
3 2015-01-01 03:00:00 983.79 983.79
4 2015-01-01 04:00:00 998.67 998.67
5 2015-01-01 05:00:00 1045.17 1045.17
6 2015-01-01 06:00:00 1085.06 1085.06
7 2015-01-01 07:00:00 1142.36 1142.36
8 2015-01-01 08:00:00 1206.03 1206.03
9 2015-01-01 09:00:00 1257.98 1257.98
total_gen_capacity_MW total_gen_res_MW total_gen_costs_€ \
0 3844.767366 1069.151457 58533.364265
1 3844.767366 1069.151457 58533.364265
2 3844.767366 1029.096172 58533.364265
3 3844.767366 1008.425544 58533.364265
4 3844.767366 1023.071677 58533.364265
5 3844.767366 1068.856000 58533.364265
6 3844.767366 1108.150370 58533.364265
7 3844.767366 1164.624394 58533.364265
8 3844.767366 1227.417593 58533.364265
9 3844.767366 1278.684192 58533.364265
max_line_loading_% max_loading_index
0 65.368569 6459
1 65.368569 6459
2 67.046971 6459
3 67.913182 6459
4 67.299425 6459
5 65.380948 6459
6 63.734596 6459
7 61.368737 6459
8 58.738544 6459
9 56.591464 6459
ac_pf converged successfully.
7.1.3 DC optimal power flow calculation with fallback mechanisms#
# Run DC optimal power flow
if not converged:
converged = False
print("AC power flow did not converge!")
print("Trying DC OPF...")
# **Step 2: First DC OPF attempt**
converged = handle_pf_results(
"dc_opf", net, time_series_short, processed_regions, pp
)
if converged:
print("Pf results on first try")
else:
# Modify generator parameters for another attempt
net.gen.loc[
(net.gen["controllable"]) & (~net.gen["type"].isin(["solar", "wind"])),
"max_p_mw",
] = net.gen.loc[
(net.gen["controllable"]) & (~net.gen["type"].isin(["solar", "wind"])),
"p_mw",
]
net.gen.loc[
(net.gen["controllable"]) & (~net.gen["type"].isin(["solar", "wind"])),
"min_p_mw",
] = 0
# **Step 3: Second DC OPF attempt**
converged = handle_pf_results(
"dc_opf", net, time_series_short, processed_regions, pp
)
if converged:
print(
"DC OPF converged setting conventional generators to minimal production"
)
else:
# Set all controllable generators to minimum production
net.gen.loc[net.gen["controllable"], "max_p_mw"] = net.gen.loc[
net.gen["controllable"], "p_mw"
]
net.gen.loc[net.gen["controllable"], "min_p_mw"] = 0
# **Step 4: Third DC OPF attempt (adjust loads)**
converged = handle_pf_results(
"dc_opf", net, time_series_short, processed_regions, pp, True
)
if converged:
print("DC OPF converged setting all generators to minimal production")
else:
net.load["min_p_mw"] = net.load["p_mw"] * 0.1 # Allow load curtailment
net.load["max_p_mw"] = net.load["p_mw"] * 1 # Keep max load unchanged
net.load["controllable"] = True # Enable load control
# **Step 5: Final attempt before failure**
converged = handle_pf_results(
"dc_opf", net, time_series_short, processed_regions, pp, True
)
if converged:
print("DC OPF converged setting variable loads")
else:
print("DC OPF did not converge!")
simple_plotly(net) # Generate diagnostic visualization
diagnose = pp.diagnostic(net) # Run diagnostics
print("Diagnostic result:", diagnose)
AC power flow did not converge!
Trying DC OPF...
Error in dc_opf: Optimal Power Flow did not converge!
time_step total_load_ref_MW total_load_res_MW \
0 2015-01-01 00:00:00 1045.47 1045.47
1 2015-01-01 01:00:00 1045.47 1045.47
2 2015-01-01 02:00:00 1004.79 1004.79
3 2015-01-01 03:00:00 983.79 983.79
4 2015-01-01 04:00:00 998.67 998.67
5 2015-01-01 05:00:00 1045.17 1045.17
6 2015-01-01 06:00:00 1085.06 1085.06
7 2015-01-01 07:00:00 1142.36 1142.36
8 2015-01-01 08:00:00 1206.03 1206.03
9 2015-01-01 09:00:00 1257.98 1257.98
total_gen_capacity_MW total_gen_res_MW total_gen_costs_€ \
0 3844.767366 1045.47 10628.404771
1 3844.767366 1045.47 10628.404771
2 3844.767366 1004.79 10058.344900
3 3844.767366 983.79 9764.066205
4 3844.767366 998.67 9972.583680
5 3844.767366 1045.17 10624.200789
6 3844.767366 1085.06 11183.190176
7 3844.767366 1142.36 11986.150614
8 3844.767366 1206.03 12878.375588
9 3844.767366 1257.98 13606.365023
max_line_loading_% max_loading_index
0 15.774209 5789
1 15.774209 5789
2 15.015744 5789
3 14.624206 5789
4 14.901638 5789
5 15.768615 5789
6 16.512351 5789
7 17.580690 5789
8 18.767796 5789
9 19.736387 5789
dc_opf converged successfully.
DC OPF converged setting conventional generators to minimal production
7.1.4 Calculating PF with scaling factor#
def calculate_scaling_factor(net):
total_gen_power = net.gen.loc[net.gen.in_service, "p_mw"].sum()
total_load_power = net.load.loc[net.load.in_service, "p_mw"].sum()
if total_gen_power > 0:
scaling_factor = total_load_power / total_gen_power
else:
scaling_factor = 1
print(f"Total generator power: {total_gen_power:.2f} MW")
print(f"Total load power: {total_load_power:.2f} MW")
print(f"Scaling factor: {scaling_factor:.3f}")
return scaling_factor
if not converged:
try:
# Calculate scaling factor
scaling_factor = calculate_scaling_factor(net)
net.gen["scaling"] = scaling_factor
# Ensure the function call includes required arguments
pf_results = time_series_pf_results(
"ac_pf", time_series_short, net, processed_regions, pp, True
)
# Check if all values in pf_results are NaN
if pf_results.isna().all().all():
print("All values in pf_results are NaN, skipping plot.")
else:
print(pf_results)
pf_res_plotly(net)
except KeyError as e:
print(f"KeyError encountered: {e}")
print("AC power flow did not converge due to a missing key!")
except ValueError as e:
print(f"ValueError encountered: {e}")
print("AC power flow did not converge due to invalid values!")
except Exception as e: # Catch-all for unexpected errors
print(f"Unexpected error: {e}")
print("AC power flow also did not converge with scaling!")
7.2 Separate Results#
# Change False to True for Separate Results - Data will be exported as excel files
separate_results = False
def run_pf_results(mode, time_series_short, processed_regions, net, pp):
"""Runs the power flow analysis and handles errors properly."""
try:
time_series_pf_results_separate_exports(
mode, time_series_short, processed_regions, net, pp
)
print(f"{mode} power flow analysis completed successfully.")
except KeyError as e:
print(f"KeyError encountered in {mode}: {e}")
except ValueError as e:
print(f"ValueError encountered in {mode}: {e}")
except Exception as e:
print(f"Unexpected error in {mode}: {e}")
print(f"{mode} power flow did not converge!")
if separate_results:
# **Execute all power flow analyses**
for mode in ["ac_pf", "dc_opf", "ac_opf"]:
run_pf_results(mode, time_series_short, processed_regions, net, pp)
8 Plotting#
8.1 Plotting Connections Between Buses and Loads / Generation#
#Plotting Processed NUTS Regions:
plot_regions(processed_regions)
# Plotting Connections Between NUTS Centroids and Buses:
plot_bus_nuts_connections(processed_regions=processed_regions, bus_geo=net.bus_geodata, net=net)
#Plotting Connections Between Generators and Buses:
plot_bus_gen_connections(network=net, bus_geo=net.bus_geodata, nuts_regions=processed_regions, plants=plants)
8.2 Plotting Results - Generation and Consumption#
#Capacities from DC OPF results - combined:
plot_power_flow_with_lines(net, res=True, sum_p_mw=True)
# Installed capacities per bus - combined
plot_power_flow_with_lines(net, res=False, sum_p_mw=True)
# Installed capacities - real locations
plot_power_flow_with_lines(net, res=False, plants=plants)